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ABSTRACT 


Ultrasonic  imaging  plays  an  important  role  as  an  adjunct  to  mammography,  with  an  emerging 
role  in  breast  cancer  screening.  Ultrasound's  real-time  nature,  lack  of  ionizing  radiation, 
and  relative  comfort  for  the  patient  make  it  an  attractive  imaging  choice.  Unfortunately, 
ultrasound  image  quality  is  often  limited. 

We  hypothesize  that  bright  scatterers  seriously  degrade  ultrasound  images  by  introducing 
image  clutter.  In  the  breast  bright  off-axis  echoes  may  originate  from  Cooper's  ligaments, 
structured  glandular  tissue,  calcification,  fat-soft  tissue  interfaces,  or  other  structures. 

While  we  initially  proposed  using  a  variant  of  the  Frost  Adaptive  Beamformer  to  reduce 
clutter,  we  have  since  discovered  that  this  technique  is  non-optimal  for  our  application. 
Extensive  literature  reviews  have  led  us  to  utilize  a  recently  proposed  method,  Spatial 
Processing  Optimized  and  Constrained  (SPOC).  In  initial  simulations  this  method  not  only 
dramatically  reduces  image  clutter,  but  also  yields  super-resolution.  We  are  actively 
.refining  this  method  while  developing  the  experimental  tools  needed  for  in  vivo  testing. 
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Introduction: 

Ultrasonic  imaging  currently  plays  an  important  role  as  an  adjunct  to  mammography  [1,  2],  Ultrasound’s  real¬ 
time  nature,  lack  of  ionizing  radiation,  and  relative  comfort  for  the  patient  make  it  an  attractive  choice  for 
applications  which  include  the  differentiation  of  fluid  filled  cysts  and  solid  masses,  differentiation  of  benign  and 
malignant  lesions,  and  guidance  of  needle  and  core  biopsy  procedures.  Recent  studies  have  even  shown  the 
potential  of  ultrasound  as  a  screening  tool,  especially  for  pre-menopausal  women  whose  radio-dense  breast 
tissue  seriously  limits  x-ray  mammography  [3],  In  both  differential  diagnosis  and  screening  applications 
however,  ultrasound  image  quality  is  limited,  with  high  levels  of  background  clutter  representing  a  significant 
problem  in  many  patients. 

While  the  cause  of  high  background  clutter  and  poor  breast  image  quality  has  not  been  determined  with 
certainty,  it  is  widely  held  that  acoustic  velocity  inhomogeneities  in  breast  tissue  cause  defocusing  of  the 
acoustic  beam.  This  distortion  manifests  itself  through  mainlobe  broadening  and  increasing  sidelobe  levels. 
Numerous  researchers,  including  the  Principal  Investigator,  have  suggested  that  this  problem,  known  as  phase 
aberration,  might  be  corrected  through  the  application  of  compensating  time  delays  [4-6],  a  combination  of 
delay  and  amplitude  corrections  [7],  or  other  more  complex  techniques  [8-10],  While  proposed  phase  aberration 
correction  methods  have  been  shown  to  have  great  potential  in  ex  vivo  or  other  non  real-time  environments, 
there  has  been  limited  evidence  of  significant  clinical  image  improvement.  The  development  of  a  real-time 
phase  aberration  correction  system  at  the  GE  Global  Research  Center  has  shown  that  real-time  phase  correction 
is  possible,  however  in  vivo  results  using  1.5-D  arrays  show  contrast  improvements  of  only  about  3  dB  in  the 
abdomen  [11,  12].  This  unimpressive  outcome  may  result  from  imperfect  algorithm  optimization,  or  perhaps  a 
lower  level  of  in  vivo  phase  aberration  than  previously  suspected.  This  latter  hypothesis  is  supported  by  recent 
phase  aberration  measurements  performed  at  Duke  University  which  indicate  in  vivo  phase  aberrations  of  only 
~25ns  RMS  (Root  Mean  Squared)  with  a  3.5  mm  FWHM  (Full  Width  at  Half  Maximum)  [13].  The  limited 
improvements  of  real-time  phase  correction,  coupled  with  low  measured  aberrations,  suggest  that  phase 
aberration  may  not  represent  the  major  source  of  breast  image  degradation. 

If  phase  aberration  is  not  the  primary  factor  limiting  breast  image  quality,  then  what  is?  We  hypothesize  that 
localized  bright  scatterers  seriously  degrade  ultrasound  images  by  introducing  broad  image  clutter.  Figure  1 
shows  single  channel  Radio  Frequency  (RF)  echo  data  obtained  from  calcifications  in  the  thyroid  of  a  human 
subject  at  Duke  University.  A  focused  transmit  beam  was  used  and  RF  data  was  acquired  from  each  element  in 
a  1.5-D  array  consisting  of  approximately  1000  elements.  Figure  1  shows  data  from  one  row  of  this  array  after 
application  of  geometrically  determined  focal  delays.  At  least  three  clear  waveforms  are  visible  in  this  data  set, 
with  each  probably  resulting  from  a  single  calcification.  Although  summation  across  channels  to  form  an  RF 
image  line  would  amplify  the  echo  coming  from  directly  in  front  of  the  array,  it  would  not  entirely  eliminate  the 
two  other  visible  targets.  These  non-focal  targets  would  appear  in  this  image  line  as  clutter,  reducing  image 
contrast.  In  addition  to  the  three  dominant  calcification  waveforms,  the  data  set  also  includes  echoes  from 
background  speckle.  These  background  echoes  also  include  discernable  off-axis  scatterers  that  undoubtedly 
generate  further  clutter  in  the  image.  Note  that  the  thyroid  data  presented  in  figure  1  is  similar  in  appearance  to 
breast  data  obtained  at  Duke.  In  the  breast  bright  off-axis  echoes  may  originate  from  Cooper’s  ligaments,  highly 
structured  glandular  tissue,  calcification,  fat-soft  tissue  interfaces,  or  other  tissue  structures. 

The  presence  of  bright  off-axis  scatterers  and  the  image  degradation  that  they  cause  is  not  surprising.  It  is  well 
known  that  the  acoustic  reflectivity  of  targets  within  the  body  covers  many  orders  of  magnitude.  It  is  precisely 
for  this  reason  that  manufacturers  employ  aggressive  apodization  to  reduce  sidelobe  levels  in  diagnostic 
ultrasound.  It  has  also  been  argued  that  harmonic  imaging  is  effective  at  improving  image  quality  because  it 
further  reduces  sidelobe  levels  and  therefore  reduces  the  spatial  spread  of  bright  targets.  The  detrimental  impact 
of  bright  scatterers  on  ultrasound  image  quality  is  recognized  in  experimental  data,  physical  intuition,  and  years 
of  experience  in  ultrasound  system  design. 
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Figure  1:  Single  channel  data  obtained  in  vivo  from  a  human  thyroid.  At  least  three  bright 
scatterers  (likely  microcalcifications)  are  visible.  While  one  target  lies  near  the  focus,  the  other 
two  are  off-axis  and  will  contribute  clutter  to  this  image  line.  Off-axis  targets  are  also  visible 
within  the  speckle  generating  background  echoes.  (Data  courtesy  Gregg  Trahey,  Duke 
University.) 

The  impact  of  a  few  bright  targets  in  an  otherwise  dim  image  has  been  well  studied  in  both  RADAR  and 
SONAR.  In  SONAR  the  detection  of  an  intentionally  stealthy  submarine  among  many  noisy  ships  requires 
separating  out  a  signal  that  is  orders  of  magnitude  below  the  signals  from  other  nearby  targets.  A  broad  variety 
of  adaptive  beamforming  algorithms  have  been  developed  for  this  scenario.  It  is  the  goal  of  this  research  project 
to  evaluate  the  potential  applicability  of  these  methods  in  medical  ultrasound. 
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Body: 

CAB  Progress: 

In  our  original  proposal  we  described  a  research  plan  to  explore  the  potential  of  one  particularly  elegant  and 
well  studied  method  for  adaptive  beamforming.  This  technique,  developed  by  Frost  in  the  early  70s,  represents 
the  beamforming  process  using  a  linear  algebra  formulation.  A  set  of  beamformer  weights,  corresponding  to  an 
individual  FIR  filter  for  each  element  of  the  array,  is  determined  adaptively  from  the  data  acquired  by  the  sensor 
array.  The  weights  are  determined  so  as  to  minimize  the  energy  in  the  beamsum  signal,  subject  to  the  constraint 
that  the  beamformer  must  exhibit  a  given  impulse  response  for  signals  received  from  the  array  look  direction. 
Over  decades  of  analysis,  the  Frost  Adaptive  Beamformer  (FAB)  has  proven  to  be  a  valuable  tool  for 
identifying  dim  targets  within  clutter  backgrounds. 

In  work  performed  prior  to  the  initiation  of  this  grant  we  built  upon  FAB  and  modified  it  to  account  for  the 
unique  characteristics  of  medical  ultrasound.  The  resultant  algorithm,  referred  to  as  the  Constrained  Adaptive 
Beamformer  (CAB),  works  in  transmit-receive  systems,  can  be  applied  to  near-field  data,  and  can  handle  the 
broad  bandwidths  and  limited  stationarity  of  ultrasound  data  [14,  15].  As  our  initial  proposal  indicated,  the  CAB 
showed  great  potential  with  both  simulation  and  experimental  ultrasound  data. 

Unfortunately  CAB  presented  significant  difficulties.  First,  because  our  implementation  of  the  CAB  required 
the  inversion  of  large  matrices  for  each  pixel  in  the  output  image,  the  computational  cost  of  implementing  CAB 
was  enormous.  Computation  of  modest  images  only  a  few  millimeters  on  a  side  took  a  day  or  more  on  a  high- 
end  PC.  While  irritating,  this  limitation  was  not  considered  fundamental.  More  serious  was  the  problem  of 
target  dependency  that  we  found  as  we  applied  CAB  to  various  data  sets.  When  applied  to  point  target  data  we 
found  that  straight-forward  implementation  of  the  algorithm  termed  Single-Iteration  CAB  (SI-CAB)  yielded 
excellent  results;  the  main-lobe  was  narrowed  and  sidelobes  were  reduced,  improving  both  spatial  resolution 
and  image  contrast.  When  we  applied  the  same  algorithm  to  data  from  a  point  target  embedded  within  a  speckle 
target  however  the  results  were  much  more  complicated.  While  the  center  of  the  point  target  was  better  focused, 
strange  behavior  was  exhibited  in  the  “tails”  or  side-lobes  of  the  point  response.  In  these  areas  the  algorithm  not 
only  eliminated  the  clutter  contribution  of  the  bright  target,  but  it  also  eliminated  the  background  speckle  signal. 
Thus  the  output  image  was  basically  black  in  the  tail  regions  of  point  target,  regardless  of  what  sort  of  target 
was  actually  present.  This  was  clearly  an  unsatisfactory  result.  In  an  attempt  to  compensate  for  this  problem  a 
previous  student  of  the  Principal  Investigator  (Jake  Mann)  developed  an  alternate  version  of  the  algorithm  that 
he  called  Multi-Iteration  CAB  or  MI-CAB.  This  algorithm  repeatedly  performed  the  SI-CAB  algorithm  on 
windows  of  data  with  partial  overlap.  It  then  averaged  output  results  in  areas  of  overlap.  Although  this  approach 
eliminated  some  of  the  artifacts  present  in  speckle  regions  it  also  showed  much  worse  performance  than  SI- 
CAB  for  point  targets.  Thus  we  were  left  with  the  highly  unsatisfactory  situation  of  requiring  different 
algorithms  for  different  target  types. 

Upon  initiation  of  this  grant  I  decided  that  it  was  critical  that  we  not  only  understand  CAB  more  thoroughly  and 
fundamentally  than  we  had  before,  but  that  we  also  explore  the  broad  array  of  adaptive  beamforming  methods 
that  had  been  developed  since  Frost’s  seminal  work  in  1972.  These  two  goals  were  far  from  independent;  the 
insights  gained  in  the  literature  review  gave  us  many  new  insights  into  the  potential  and  limitations  of  CAB. 

In  our  review  of  the  literature  we  came  to  realize  more  deeply  that  the  FAB  was  designed  for  problems  that 
were  fundamentally  different  than  that  we  face  in  ultrasound  imaging.  We  were  reminded  that  the  FAB  was 
designed  to  address  situations  where  the  targets  were  in  the  far-field  of  the  imaging  array  and  where  the  acoustic 
environment  was  effectively  narrowband.  These  assumptions  we  had  been  aware  of  and  had  addressed  carefully 
in  the  design  of  CAB.  What  we  had  been  less  clear  about  was  the  assumption  made  by  FAB  that  there  would  be 
multiple,  statistically  independent,  signal  realizations  available  from  any  given  field  of  targets.  In  medical 
ultrasound  this  assumption  is  completely  wrong.  We  can  only  get  one  realization  of  the  signals  from  any 
particular  target  arrangement.  (This  is  not  strictly  true  if  compounding  is  applied,  but  in  normal  imaging  this  is  a 
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fair  assessment.)  This  difference  proves  to  be  huge  and  is  the  main  reason  that  we  found  such  poor  behavior 
when  applying  CAB  to  speckle  targets.  While  compounding  could  theoretically  be  used  to  get  multiple 
independent  looks  at  the  target,  a  careful  review  of  the  literature  indicated  that  practical  compounding  schemes 
could  never  yield  enough  independent  looks  to  provide  statistically  robust  performance.  Thus  a  careful  literature 
review  followed  by  a  series  of  well  designed  computer  simulations  showed  us  that  CAB  was  a  fundamentally 
flawed  algorithm  for  medical  ultrasound  imaging.  Clearly  we  would  need  to  identify  another  algorithmic 
foundation  for  this  research  effort. 

Review  of  the  Adaptive  Beamforming  Literature: 

Our  review  of  the  adaptive  beamforming  literature  covered  dozens  of  papers  and  for  many  of  these  we  followed 
up  by  coding  the  described  algorithms  in  MATLAB  and  testing  their  performance  on  simulation  ultrasound 
data.  In  the  process  of  reviewing  this  area  we  found  it  particularly  helpful  to  categorize  algorithms  depending 
upon  their  assumptions  and  limitations.  There  were  obvious  differences  between  algorithms  that  were 
developed  for  near-field  application  and  those  intended  for  far-field  implementation.  Likewise,  algorithms 
intended  for  narrow-band  signals  could  be  readily  differentiated  from  those  intended  for  broad-band  signals.  For 
both  of  these  differentiators  we  found,  as  we  had  with  the  development  of  CAB,  that  it  was  often  (although  not 
always)  straightforward  to  identify  algorithm  modifications  that  would  move  them  from  one  domain  to  another. 
A  more  fundamental  difference,  that  was  much  less  clear  at  the  beginning  of  our  work,  was  the  importance  of 
the  number  of  available  realizations  of  the  data.  Specifically,  some  algorithms  assumed  that  data  from  the 
targets  could  be  acquired  repeatedly  with  each  acquisition  offering  a  new  statistical  realization  of  the  system. 
Other  algorithms  assumed  that  only  one  look  at  the  data  was  available.  This  difference  proves  to  be  critical  in 
our  problem. 

To  summarize,  the  problem  of  adaptive  beamforming  in  breast  imaging  can  generally  be  stated  as  utilizing 
broad-band  signals,  operating  in  the  near-field  of  the  transducer  arrays,  and  operating  on  only  one  look  at  the 
target. 

In  the  following  sections  we  summarize  our  findings  in  the  adaptive  beamforming  literature  and  classify 
algorithms  in  terms  of  their  assumptions  and  potential  for  our  problem. 

Statistical  beamformers 

The  following  beamformers  make  use  of  the  second  order  statistic  of  the  data  in  order  to  determine  filter 
weights  for  the  rejection  of  off-axis  signals.  They  have  been  well  characterized  and  some  have  been  known  for 
decades. 

Capon  Beamformer  (also  know  as  the  Minimum  Variance  Distortionless  Response,  or  MVDR, 
beamformer) 

This  algorithm  solves  for  a  weight  vector  that  minimizes  the  energy  of  the  beamformer’ s  output  given  a 
linear  constraint.  The  sum  of  the  weights  is  equal  to  some  number  g  (generally,  g=l). 

The  Capon  BF  requires  narrowband  signals. 

The  main  drawback  of  this  algorithm  is  that  it  requires  knowledge  of  the  second  order  statistics  of  the 
data  (i.e.,  autocorrelation  matrix).  The  autocorrelation  matrix  cannot  be  reliably  determined  without 
acquiring  many  statistically  independent  looks  at  the  target.  This  is  not  possible  in  medical  ultrasound, 
because  even  compounding  does  not  yield  enough  independent  looks  to  produce  a  useful  image. 

J.  Capon,  “High  resolution  frequency-wavenumber  spectrum  analysis,”  Proc.  IEEE,  vol.  57,  pp.  1408- 
1418,  1969. 
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Frost  Beamformer 

Also  in  this  case,  a  weight  vector  is  chosen  to  minimize  the  output  energy  given  a  linear  constraint  on 
the  weights.  The  main  difference  with  the  Capon  BF  is  that  this  algorithm  uses  K  taps  in  time,  therefore 
allowing  the  use  of  broadband  signals. 

Again,  the  autocorrelation  matrix  of  the  data  is  required.  Furthermore,  significant  problems  occur  if  the 
signal  of  interest  and  noise  (off-axis  signal)  are  correlated  (i.e.,  they  have  the  same  spectral  response). 
This  is  a  common  occurrence  in  medical  ultrasound. 

O.  L.  Frost  III,  “An  Algorithm  for  Linearly  Constrained  Adaptive  Array  Processing,”  Proc.  IEEE ,  vol. 

60,  no.  8,  pp.  926-935,  1972. 

Generalized  Sidelobe  Cancellor  (Griffith’s  formulation) 

Equivalent  to  the  Frost  algorithm  in  its  formulation  and  assumptions.  The  main  advantage  of  the  GSC  is 
that  it  transforms  the  minimization  from  constrained  to  unconstrained,  therefore  increasing  the  degree  of 
freedoms  of  the  beamformer.  This  is  accomplished  using  a  subtractive  process  on  the  array  elements. 
The  GSC  has  many  of  the  same  limitations  as  the  Frost  BF,  but  is  somewhat  easier  to  implement. 

L.  J.  Griffith  and  C.  W.  Jim,  “An  Alternative  Approach  to  Linearly  Constrained  Adaptive 
Beamforming,”  IEEE  Trans.  Antennas  Propagat.,  vol.  AP-30,  no.  1 ,  pp.  27-34,  1982. 

Duvall  Beamformer 

Equivalent  to  the  Frost  beamformer  in  its  formulation  and  assumptions.  This  algorithm,  however,  is 
designed  to  get  around  the  problem  of  correlated  signal  and  noise.  This  is  critical  for  medical  ultrasound, 
where  clutter  usually  results  from  reflection  of  the  same  pulse  as  that  used  to  image  the  target.  The 
Duvall  BF  uses  two  BFs;  the  first  one  is  a  slave  BF,  whereas  the  second  is  connected  to  the  array 
elements  through  a  subtractive  process.  The  weights  calculated  in  the  latter  BF  are  then  copied  in  the 
slave  BF  to  generate  an  output.  While  this  algorithm  is  in  some  ways  attractive,  it  still  requires  enough 
target  realizations  to  generate  a  good  autocorrelation  matrix;  a  requirement  that  cannot  be  met  in 
medical  ultrasound. 

B.  Widrow,  K.  M.  Duvall,  R.  P.  Gooch,  and  W.  C.  Newman,  “Signal  Cancellation  Phenomena  in 
Adaptive  Antennas:  Causes  and  Cures,”  IEEE  Trans.  Antennas  Propagat.,  vol.  AP-30,  no.  3,  pp.  469- 
478, 1982. 

Reduced  Rank  Beamformers 

The  basic  idea  here  is  to  save  computation  time  by  calculating  a  reduced  rank  autocorrelation  matrix.  All 
considerations  made  for  the  Capon  BF  still  apply  in  this  case.  Principal  Component  (PC),  and  Multi- 
Stage  Wiener  Filter  (MSWF)  are  example  of  such  algorithms. 

W.  F.  Gabriel,  “Using  Spectral  Estimation  Techniques  in  Adaptive  Processing  Antenna  Systems,”  IEEE 
Trans.  Antenna  Propagat.,  vol.  AP-34,  no.  3,  pp.  291-300,  1986. 

S.  M.  Kogon,  “Robust  adaptive  beamforming  for  passive  sonar  using  eigenvector/beam  association  and 
excision,”  Sensor  Array  and  Multichannel  Signal  Processing  Workshop  Proceedings,  pp.  33-37,  2002. 

J.  S.  Goldstein,  I.  S.  Reed,  and  L.  L.  Scharf,  “A  Multistage  Representation  of  the  Wiener  Filter  Based  on 
Orthogonal  Projections,”  IEEE  Trans.  Inform.  Theory,  vol.  44,  no.  7,  pp.  2943-2959,  1998. 

M.  L.  Honig  and  J.  S.  Goldstein,  “Adaptive  Reduced-Rank  Interference  Suppression  Based  on  the 
Multistage  Wiener  Filter,”  IEEE  Trans.  Commun.,  vol.  50,  no.  6,  pp.  986-994,  2002. 
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Deterministic  Beamformers 

The  following  beamformers  do  not  make  use  of  any  statistics,  therefore  are  suitable  to  use  with  only  a  single¬ 
snapshot  of  data.  In  general,  this  class  of  beamformer  is  much  better  suited  to  the  problems  we  face  in 
ultrasound  imaging.  These  algorithms  are  relatively  new  and  research  is  progressing  rapidly  in  this  domain. 

Adaptive  Single  Snapshot  Beamformer 

This  algorithm  subdivides  the  array  into  G  groups  of  K  elements  (the  groups  can  overlap),  and  rearranges 
the  data  into  a  GxK  matrix.  For  every  steering  direction,  this  beamforming  scheme  searches  for  the  scalar  a 
that  reduced  the  rank  of  the  data  matrix.  This  scalar  represents  the  magnitude  of  the  signal  coming  from  that 
particular  look  direction.  This  algorithm  is  based  on  the  following  observation:  each  source  in  the  far  field 
contribuites  as  1  in  the  rank  of  the  GxK  matrix.  In  other  words,  if  we  assume  L  far  field  sources,  the  rank  of 
the  data  matrix  will  be  equal  to  L. 

The  ASS  beamformer  is  designed  for  narrowband,  far  field  signals.  However,  it  can  also  be  implemented  in 
the  frequency  domain  for  broadband  data.  It  is  apparently  not  flexible  enough  to  be  applied  to  near-field 
signals  however. 

M.  E.  Ali  and  F.  Schreib,  “Adaptive  Single  Snapshot  Beamforming:  A  New  Concept  for  the  Rejection  of 
Nonstationary  and  Coherent  Interferes,”  IEEE  Trans,  Signal  Processing,  vol.  40,  no.  12,  pp.  3055-3058, 
1992. 

Deterministic  Null  Steering  (DNS) 

The  main  idea  behind  this  algorithm  is  to  project  the  data  into  a  sub-space  that  is  orthogonal  to  the  space 
spanned  by  the  off-axis  signals.  In  order  to  do  that,  the  algorithms  required  knowledge  of  the  DOAs  of  the 
signals  that  need  to  be  null  out.  (If  DOA’s  have  to  be  estimated  from  the  data,  this  algorithm  should  be 
considered  a  reduced  rank  adaptive  BF). 

In  the  narrowband  case  the  algorithm  is  easily  implemented,  since  the  DOAs  are  one-dimensional  steering 
vectors.  However,  in  the  broadband  case,  the  computational  complexity  increases  significantly  since  we 
have  to  deal  with  FIR  filters  (preferably  matched  filters)  for  every  channel. 

An  advantage  of  this  algorithm  is  that  it  requires  no  far-field  assumption.  Unfortunately  algorithm 
performance  is  highly  sensitive  to  mismatches  in  DOA.  In  simulations  we  concluded  that  this  sensitivity 
was  so  great  that  this  algorithm  was  unlikely  to  perform  well  with  experimental  ultrasound  data. 

H.  Subbaram  and  K.  Abend,  “Interference  Suppression  Via  Orthogonal  Projections:  A  Performance 
Analysis,”  IEEE  Trans.  Antennas  and  Propag.,  vol.  41,  no.  9,  pp.  1187-1 194,  1993. 

H.  L.  Van  Trees,  “Detection,  Estimation,  and  Modulation  Theory;  Part  IV:  Optimum  Array  Processing,” 
John  Wiley  &  Sons,  New  York,  2002. 

Source  Localization  Algorithms 

This  class  of  algorithms  takes  a  completely  different  approach.  A  signal  model  is  first  generated  of  the  field 
produced  by  a  series  of  hypothetical  sources.  Through  some  recursive  algorithms,  the  sensors’  output  is 
matched  to  the  signal  model  to  solve  for  the  position  and  intensity  of  the  real  sources.  Statistics  are  not 
required  and  both  near-field/far-field  and  broad/narrow  band  cases  can  be  modeled  (for  the  broadband  case 
the  system  is  usually  modeled  in  the  frequency  domain). 

The  main  problems  with  this  class  of  algorithm  are  the  high  computational  cost,  the  derivation  of  an 
accurate  signal  model,  and  possibly  convergence.  These  limitations  not  withstanding,  we  have  decided  to 
focus  our  efforts  on  this  class  of  algorithms. 

R.  Bethel,  B.  Shapo,  and  H.  L.  Van  Trees,  “Single  Snapshot  Spatial  Processing:  Optimized  and 
Constrained,”  Sensor  Array  and  Multichannel  Signal  Processing  Workshop  Proceedings,  pp.  508-512, 

2002. 
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D.  M.  Malioutov,  M.  Cetin,  J.  W.  Fisher  III,  and  A.  S.  Willsky,  “Superresolution  Source  Localization 
Through  Data- Adaptive  Regularization,”  Sensor  Array  and  Multichannel  Signal  Processing  Workshop 
Proceedings,  pp.  194-198,  2002. 

After  performing  extensive  literature  searches  and  implementing  and  testing  a  number  of  the  algorithms 
described  above,  we  concluded  that  the  algorithm  described  in  “Single  Snapshot  Spatial  Processing:  Optimized 
and  Constrained”  offered  the  best  possibility  for  success  given  the  parameters  of  in  vivo  ultrasound  imaging 

SPOC  Progress: 

Most  of  the  Adaptive  Beamforming  algorithms  previously  described  fail  when  applied  to  medical 
ultrasound  data.  This  can  be  attributed  to  some  or  all  of  the  following  factors:  we  are  operating  in  a  near-field 
scenario,  our  signals  are  broadband,  and  we  have  limited  statistical  information  available. 

SPOC  was  first  developed  by  Van  Trees  et  al..  A  similar  approach  has  also  been  discussed  by  Malioutov 
et  al..  This  algorithm,  briefly  described  below,  can  be  successfully  applied  to  medical  ultrasound  data. 

The  algorithm: 

•  The  region  of  interest  (ROI)  is  subdivided  in  a  grid  of  hypothetical  sources  at  arbitrary  positions,  as 
shown  in  figure  2.  Finer  grid  sampling  yields  finer  final  resolution. 

•  For  each  point  in  the  ROI,  we  calculate  the  hypothetical  field  received  by  the  array  for  that  specific 
point.  These  responses  are  included  into  the  array  manifold  matrix  V  of  dimensions  LxP,  where  L  and  P 
are  the  numbers  of  hypothetical  sources  in  the  range  and  lateral  dimensions,  respectively. 

•  The  data  vector  x  =  [x,  x2  ...  xN]  is  received  by  the  N-element  linear  array.  This  vector  can  be  modeled 
as  x  =  V*f.  In  this  expression,  f  is  the  LPxl  signal  vector,  whose  elements  are  the  amplitudes  of  sources 
located  in  the  ROI. 

•  Given  x  and  V,  SPOC  recursively  estimates  the  signal  vector  f. 


Figure  2.  Schematic  diagram  of  SPOC. 


A  series  of  computer  simulations  were  performed  in  Matlab  to  test  the  potential  of  SPOC. 

Point  Target. 

We  simulated  a  33  element  linear  array  operating  at  5  MHz  with  element  spacing  of  150  pm.  A  point  target  was 
placed  in  front  of  the  transducer  at  a  depth  of  20. 1  mm. 

Hypothetical  sources  were  placed  every  20_m  in  range  and  every  100pm  in  azimuth. 
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Figure  3.  Results  obtained  simulating  a  point  target.  The  point  target  is  depicted  in  (a);  (b)  shows  conventional 
delay-and-sum  beamforming,  whereas  images  (c)  to  (f)  show  SPOC  with  different  levels  of  electronic  noise. 


Wire  Targets. 

We  simulated  a  33  element  linear  array  operating  at  5  MHz  with  element  spacing  of  150  pirn.  A  series  of  wires 
were  placed  in  front  of  the  transducer  at  various  depths  and  lateral  positions. 

Hypothetical  sources  were  placed  every  20/*  m  in  range  and  every  100/*m  in  azimuth. 
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Figure  4.  Results  obtained  simulating  a  series  of  wire  targets.  The  wires  are  depicted  in  (a);  (b)  shows 
conventional  delay-and-sum  beamforming,  whereas  (c)  show  SPOC. 
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Anechoic  Cyst. 

We  simulated  a  33  element  linear  array  operating  at  5  MHz  with  element  spacing  of  150  //m.  An  anechoic  cyst 
in  a  scattering  environment  was  placed  in  front  of  the  transducer. 

Hypothetical  sources  were  placed  every  20// m  in  range  and  every  100/im  in  azimuth. 
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Figure  5.  Results  obtained  simulating  an  anechoic  cyst  surrounded  by  scatterers.  The  cysts  is  depicted  in  (a);  (b) 
shows  conventional  delay-and-sum  beamforming,  whereas  (c)  show  SPOC. 

Simulation  Tool  Development: 

The  performance  of  the  SPOC  algorithm  depends  greatly  upon  the  quality  of  the  system  model  it  applies.  In  our 
initial  simulations  (presented  above)  we  modeled  the  system  using  analytical  methods,  or  using  the  well  known 
FIELD  program  written  by  Jprgen  Jensen.  Unfortunately  FIELD  is  quite  slow  when  computing  full  3D  spatial 
impulse  responses.  Further,  because  FIELD  works  entirely  on  sampled  data  sets  it  is  prone  to  artifacts  from 
undersampling.  We  have  implemented  the  Tupholme-Stepanishen  method  (the  core  approach  used  in  FIELD)  in 
a  new  piece  of  code  we  call  DELFI.  The  DELFI  code  uses  cubic  spline  representations  of  the  transmitted  pulse, 
and  the  transmit  and  receive  spatial  impulse  responses.  This  approach  avoids  the  potential  sampling  difficulties 
of  FIELD.  It  is  also  significantly  faster  at  computing  space-space-space  responses  at  an  instant  in  time.  These 
sort  of  responses  are  critical  in  much  of  our  research  and  the  25  fold  increase  in  speed  for  DELFI  is  of  great 
significance.  We  have  submitted  a  paper  describing  the  algorithm  and  DELFI  code  [16].  Once  this  paper  is 
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accepted  for  publication  we  will  place  the  1000  lines  of  DELFI  source  code  on  the  Mathworks  public  access 
web  site. 

A  limitation  of  the  DELFI  code  is  that  it  utilizes  the  same  far-field  approximation  used  by  FIELD.  This 
assumption  requires  intensive  subsampling  of  the  transducer  array  when  the  point  of  interest  lies  in  the  near¬ 
field  of  the  transducer.  We  are  currently  refining  the  algorithm  to  be  more  accurate  in  the  near-field.  The 
modifications  currently  being  explored  would  not  have  a  significant  impact  on  computational  complexity.  We 
anticipate  completing  and  testing  this  refined  algorithm  in  the  coming  year.  We  will  submit  a  publication 
describing  this  work  as  it  develops. 

Experimental  Platform  Development: 

Both  the  CAB  algorithm  described  in  our  original  proposal  and  the  SPOC  algorithm  which  we  are  now 
exploring  require  access  to  individual  channel  data  within  the  ultrasound  beamformer.  Acquisition  of  such  data, 
especially  in  real-time,  is  a  daunting  task.  We  are  currently  working  with  three  different  experimental  platforms. 
Each  of  these  has  limitations,  but  together  their  capabilities  are  quite  expansive. 

Our  first  experimental  system  is  a  clinical  GE  Logiq  700  MR.  This  scanner  acquires  summed  complex  I/Q  data 
under  control  of  a  set  of  custom  pulse  sequences  and  system  tools.  We  are  able  to  obtain  single  channel  data  by 
altering  system  apodization  on  successive  transmissions.  Acquisition  of  a  complete  single  channel  data  set 
including  transmit  focusing  is  extremely  slow  with  this  method,  however  limited  data  sets  can  be  acquired 
rapidly.  In  one  strategy  we  translate  a  single  transmit  element  and  a  single  receive  element  across  the  array.  This 
allows  synthetic  aperture  (SA)  focusing  and  CAB  to  be  performed  offline.  Surprisingly  this  approach  exhibits  a 
single  channel  signal  to  noise  ratio  (SNR)  of  approximately  36  dB.  Unfortunately  tissue  motion  during 
acquisition  can  seriously  degrade  images  formed  using  this  method.  We  have  recently  implemented  an  alternate 
acquisition  strategy  which  utilizes  a  focused  transmit  beam  and  a  moving  single  element  receive  aperture  to 
enable  single  channel  acquisition.  This  approach  will  allow  us  to  acquire  a  full  set  of  single  channel  data  for  a 
few  dozen  image  lines  in  the  same  period  of  time  that  would  be  required  to  acquire  a  normal  frame  of  image 
data.  While  still  somewhat  slow,  this  approach  should  be  less  susceptible  to  tissue  motion.  We  are  currently 
testing  this  method. 

Our  second  experimental  system  is  a  Philips  SONOS  5500.  This  system  also  operates  under  custom  control  and 
acquires  single  channel  data  through  system  apodization.  We  are  expanding  upon  this  capability  by  modifying 
this  system  to  acquire  data  in  real-time  and  in  parallel  from  all  128  beamformer  channels.  Once  complete  we 
will  be  able  to  acquire  single  channel  data  for  up  to  3.2  seconds.  Note  that  single  channel  acquisition  with  this 
system  will  be  in  real-time  and  will  include  full  use  of  transmit  focusing.  In  the  past  year  we  have  designed  and 
tested  the  final  circuit  board  needed  for  this  system  modification.  Unfortunately  testing  revealed  a  fundamental 
flaw  in  that  design.  We  have  since  redesigned  the  board  and  are  now  assembling  it  for  test.  We  think  it  likely 
that  this  system  will  become  functional  in  the  first  quarter  of  2006. 

A  third  experimental  system,  not  described  in  the  original  proposal,  is  a  fully  custom  system  developed  in  a 
collaboration  with  two  other  investigators  at  the  University  of  Virginia  (John  A.  Hossack  and  Travis  N. 
Blalock.)  A  second  generation  of  this  system  is  already  being  assembled.  This  second  generation  system  will 
utilize  a  3600  element  2D  array  and  will  operate  at  a  5.0  MHz  center  frequency.  Data  from  all  channels  will  be 
acquired  in  parallel  in  real-time,  however  only  four  real  samples  will  be  acquired  by  each  channel.  Although  the 
data  acquisition  of  this  system  is  certainly  limited  in  some  ways,  we  believe  that  this  system  will  provide  an 
excellent  testbed  for  SPOC.  In  the  coming  months  as  this  system  becomes  operational  we  will  acquire  data  to 
test  the  viability  of  SPOC  on  2D  arrays.  I 
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Key  Research  Accomplishments: 

•  Identification  of  a  fundamental  weakness  of  the  CAB  algorithm  for  ultrasound  imaging 

•  Extensive  review  of  adaptive  beamforming  literature 

•  MATLAB  implementation  and  evaluation  of  a  number  of  adaptive  beamforming  algorithms 

•  Identification  of  the  SPOC  adaptive  beamforming  algorithm 

•  MATLAB  implementation  of  SPOC 

•  Extensive  testing  of  SPOC  for  point  targets  and  speckle  generating  targets 

•  Tested  sensitivity  of  SPOC  to  electronic  noise  (not  excessively  sensitive) 

•  Tested  sensitivity  of  SPOC  to  the  effects  of  frequency  dependent  attenuation  (some  sensitivity,  can  be 
compensated  for  through  modest  algorithm  adjustment) 

•  Conceived  of  and  implemented  the  DELFI  simulation  tool 

•  Submitted  a  paper  describing  DELFI 

•  Identified  problems  with  the  SONOS  experimental  platform  under  development 

•  Revised  design  of  the  SONOS  platform 

•  Submitted  one  patent  disclosure  on  the  use  of  SPOC  in  medical  ultrasound 

•  Submitted  multiple  abstracts  to  scientific  meetings  describing  the  use  of  SPOC  in  medical  ultrasound 


Reportable  Outcomes: 

Papers: 

Walker,  W.F.,  “A  Spline  Based  Approach  for  Computing  Spatial  Impulse  Responses,”  submitted  to  IEEE  Trans.  Ultrason. 
Ferroelec.  Freq.  Contr.,  Feb.,  2005. 

Conference  Abstracts: 

Viola,  F.  and  W.F.  Walker,  “Adaptive  Signal  Processing  in  Medical  Ultrasound  Beamforming,”  accepted  to  the  2005  IEEE 
Ultrasonics  Symposium. 

Viola,  F.,  and  W.F.  Walker,  “Adaptive  Beamforming  for  Medical  Ultrasound  Imaging,”  U.S.  Dept,  of  Defense  Breast  Cancer 
Research  Program  Era  of  Hope  2005  Meeting,  June  2005. 

Patent  Disclosures: 

“Adaptive  Beamformiong  for  Medical  Ultrasound  Imaging,”  F.  Viola,  and  W.F.  Walker,  patent  disclosure  filed  2005. 

Software: 

DELFI  -  A  cubic  spline  based  code  for  simulating  spatial  impulse  responses.  Currently  over  1000  lines  of  C  code. 
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V 


Conclusions: 

In  our  first  year  of  funding  we  have  identified  a  fatal  flaw  in  the  Constrained  Adaptive  Beamforming  algorithm 
that  was  to  be  the  core  of  our  research  effort  in  this  grant.  We  have  since  explored  the  adaptive  beamforming 
literature  and  identified  an  algorithm  which  is  much  better  suited  to  our  problem:  Spatial  Processing  Optimized 
and  Constrained.  This  algorithm  is  robust  to  noise  and  other  common  difficulties  with  ultrasound  data.  In  the 
coming  year  we  will  test  this  algorithm  on  experimental  data  and  obtain  in  vivo  data  allowing  us  to  test  the 
hypothesis  that  off  axis  scattering  is  a  major  source  of  image  degradation  in  the  breast. 
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Abstract: 

This  paper  describes  an  analytical  approach  for  computing  the  two-way  far-field  spatial 
impulse  response  from  rectangular  transducer  elements  under  arbitrary  excitation.  The 
described  approach  determines  the  response  as  the  sum  of  polynomial  functions,  making 
computational  implementation  quite  straightforward.  The  proposed  algorithm  was 
implemented  as  a  C  routine  under  Matlab  and  results  were  compared  to  theory  and  to 
those  obtained  under  similar  conditions  from  the  well  established  FIELD  II  program.  All 
results  were  in  excellent  agreement.  The  proposed  method  is  quite  efficient  for 
computing  spatial  sensitivity  functions  at  a  single  instant  in  time;  an  application  for 
which  FIELD  II  is  not  well  optimized.  Under  the  specific  conditions  tested  here,  the 
proposed  algorithm  was  approximately  25  times  faster  than  FIELD  II  for  computing 
spatial  sensitivity  functions,  with  no  loss  in  quality. 

Introduction: 

The  design  of  modern  phased  array  ultrasonic  imaging  systems  relies  heavily  on 
the  use  of  computer  simulations.  This  is  necessary  because  the  broadband  and  near-field 
nature  of  most  clinical  imaging  environments  severely  limits  the  utility  of  the  Fraunhofer 
approximation  [1]  and  other  theoretical  methods.  Furthermore,  the  high  degree  of 
optimization  of  modem  systems  makes  even  small  deviations  from  theory  significant.  For 
example,  if  the  system  designer  is  concerned  with  the  array  sensitivity  pattern  down  80 
dB  from  the  main-lobe,  then  a  deviation  from  theory  of  only  0.1%  (-60dB)  will 
significantly  affect  performance  and  make  optimization  to  the  desired  level  impossible. 
Clearly,  highly  accurate  simulation  tools  are  required  to  guide  the  selection  of  transducer 
geometry,  apodization,  operating  frequency,  and  other  parameters. 

,  While  most  researchers  and  designers  would  agree  upon  the  need  for  accurate 
simulation,  the  appropriate  approach  depends  upon  the  specific  problem  of  interest  and 
the  parameters  that  are  most  significant  to  that  problem.  For  example,  in  cases  where 
details  of  transducer  vibration  and  crosstalk  are  of  interest,  a  computationally  costly  finite 
element  analysis  may  be  required  to  capture  the  most  relevant  behavior.  For  such 
problems  the  highly  optimized  PZFlex  (Weidlinger  Associates,  Inc.  New  York,  NY) 
package  is  widely  used  [2].  In  other  cases,  where  the  detailed  transducer  response  is  of 
less  interest,  but  the  propagation  medium  is  inhomogeneous  or  multiple  scattering  occurs, 
the  more  computationally  efficient  Finite  Difference  Method,  such  as  that  implemented  in 
Wave2000  (CyberLogic  Inc.,  New  York,  NY),  may  be  used. 

For  the  vast  majority  of  situations,  where  the  propagation  medium  can  be 
considered  homogeneous  or  inhomogeneities  can  be  modeled  as  a  near-field  thin  phase 
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screen,  Stepanishen’s  method  [3],  as  implemented  in  Jprgen  Jensen’s  FIELD  II  program 
[4]  has  become  a  de  facto  standard.  This  approach  determines  the  spatial  impulse 
response  of  each  transmit  element,  convolves  this  with  the  spatial  impulse  response  of 
each  receive  element,  convolves  this  result  with  the  transmitted  pulse,  and  convolves  this 
with  the  transmit  and  receive  electromechanical  impulse  responses  to  determine  the  two- 
way  temporal  response  at  a  point  in  space.  This  technique,  as  implemented  in  Jensen’s 
code,  has  been  highly  refined  over  roughly  a  decade  of  development  so  that  it  is 
extremely  efficient  and  available  in  a  compiled  form  on  a  variety  of  computer  platforms 
(http://www.es.oersted.dtu.dk/staff/jaj/field/).  By  computing  the  temporal  signals 
returned  from  various  target  locations,  FIELD  II  readily  models  common  experimental 
situations. 

While  the  temporal  response  returned  by  FIELD  II  provides  an  excellent  parallel 
to  experiment,  recent  theoretical  work  by  Zemp  [5]  and  Walker  [6]  highlights  the 
importance  of  considering  the  full  four  or  five  dimensional  system  response.  In  a 
previous  paper  (6],  this  author  discussed  a  general  signal  model  describing  the  system 
response  as  a  function  of  three  spatial  dimensions  and  time.  While  much  of  this  detail  is 
hidden  experimentally,  the  consideration  of  the  full  dimensionality  of  the  system 
response  yields  insights  and  offers  paths  for  analysis  that  are  not  apparent  in  the  more 
conventional  two  dimensional  (space,  time)  view  of  the  system.  Zemp  [5]  carried  this 
concept  one  step  further,  including  another  dimension  for  image  line  index,  thereby 
further  clarifying  system  behavior.  Our  laboratory  has  recently  applied  these  frameworks 
to  derive  a  general  resolution  metric  that  allows  quantitative  comparison  of  system 
performance,  even  when  the  individual  impulse  responses  of  those  systems  are  very 
different  in  structure  [7],  Interestingly,  this  new  resolution  metric  is  based  upon  the 
system  response  throughout  space  at  a  single  instant  in  time;  a  form  of  the  impulse 
response  that  is  not  naturally  determined  by  FIELD  II  Although  such  responses  can  be 
computed  by  sampling  the  temporal  responses  generated  by  FIELD  II,  this  approach  is 
extremely  costly  in  terms  of  both  computation  and  storage. 

In  this  paper  we  describe  a  new  approach  to  computing  spatial  impulse  responses 
that  directly  determines  the  response  throughout  space  at  a  single  instant  in  time.  This 
approach  is  complementary  to  FIELD  II,  simply  yielding  responses  in  a  different  set  of 
dimensions.  Because  results  from  this  code  are  predictive  of  system  performance  and  are 
a  permutation  of  the  data  available  from  FIELD  II,  we  name  this  code  DELFI.  In  this 
paper  we  describe  the  theoretical  underpinnings  of  the  DELFI  code,  describe 
implementation,  and  validate  the  code  through  comparisons  with  FIELD  II  and  analytical 
predictions.  We  discuss  the  relative  computational  efficiency  of  DELFI  and  discuss 
future  directions  for  development  and  refinement. 

Theory:  ; 

We  begin  our  derivation  by  considering  the  general  approach  employed  by  Jensen  in  the 
FIELD  II  program  [4].  We  consider  the  system  response  for  a  specific  transmit-receive 
element  pair  to  be  a  four  dimensional  function  of  space  and  time: 

p(x,y,z,t)  =  e(t)*ml{t)*mr  (t)*h,  (x,y,z,t)*hr(x,y,z,t)  (1) 
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where  x,  y,  and  z  are  the  three  spatial  dimensions,  t  is  the  time  for  a  given  line 
(proportional  to  range  in  the  beamformed  image),  p(x,y,z,t)  is  the  system  point  spread 
function  (psf),  <?(*)  is  the  electrical  excitation  of  the  transmit  element,  mt(t )  and  mr(t ) 
are  the  electromechanical  transfer  functions  of  the  transmit  and  receive  elements 
respectively,  h,(x,y,z,t )  and  hr(x,y,z,t)  are  the  spatial  impulse  responses  of  the 
transmit  and  receive  elements  respectively,  and  *  indicates  convolution  in  the  time 

t  < 

dimension.  In  typical  systems  the  excitation  and  the  transmit  and  receive 
electromechanical  transfer  functions  are  assumed  constant  for  all  elements  of  the  array. 
Thus  we  can  convolve  these  terms  together  before  computing  the  overall  response  with 
little  loss  in  generality.  Performing  this  step  we  simplify  (1)  to  yield: 

p(x,y,z,t)  =  emmlr(t)*hl(x,y,z,t)*hr(x,y,z,t)  (2) 

where  emmlr(t )  is  the  combined  effect  of  the  excitation  and  the  transmit  and  receive 
electromechanical  transfer  functions  and  can  be  represented  mathematically  as 
emmlr(t)=  e(t)*mt(t)*mr(t).  While  FIELD  II  computes  expressions  (1)  or  (2)  using 

sampled  versions  of  each  of  the  component  signals,  we  take  an  alternate  approach, 
instead  using  analytical  expressions  for  these  functions. 

The  utility  of  an  analytical  approach  depends  upon  the  choice  of  expressions  used: 
they  must  be  general  enough  to  include  all  relevant  cases,  but  they  must  be  constrained  in 
such  a  way  to  guarantee  the  presence  of  an  analytical  solution.  Since  (2)  allows  for 
consideration  of  most  practically  interesting  cases  and  requires  two  convolutions  (rather 
than  the  four  convolutions  of  (1))  we  build  our  algorithm  upon  this  expression.  Further 
simplification  can  be  made  by  assuming  that  the  point  of  interest  lies  in  the  far-field  of 
both  the  transmit  and  receive  elements.  This  is  not  an  onerous  assumption  since  cases 
where  the  response  would  lie  in  the  near-field  of  a  physical  element  can  be  readily 
modeled  using  a  superposition  of  computational  elements  for  which  the  far-field 
assumption  is  valid.  We  further  simplify  the  problem  by  assuming  that  the  elements  are 
rectangular. 

Following  these  assumptions,  and  once  again  drawing  upon  the  methodology  of 
Jensen  [4],  we  recognize  that  the  one-way  spatial  impulse  response  of  an  element  takes 
on  one  of  three  functions.  If  the  field  point  lies  on  a  line  perpendicular  to  the  element  face 
and  passing  through  its  center  then  the  spatial  impulse  response  as  a  function  of  time  is 
simply  a  delta  function,  as  shown  in  the  left  panel  of  figure  1 .  If  the  field  point  does  not 
fulfill  the  first  condition,  but  instead  lies  upon  one  of  two  planes  passing  through  the 
element  center  and  perpendicular  to  the  element  edges  then  the  spatial  impulse  response 
in  time  is  a  rectangle  function,  as  depicted  in  the  central  panel  of  figure  1.  Finally,  if  the 
field  point  fulfills  neither  of  the  above  conditions  then  the  spatial  impulse  response  in 
time  is  a  trapezoid  function,  as  shown  in  the  right  panel  of  figure  1.  These  possible  one¬ 
way  spatial  impulse  responses  are  summarized  mathematically  below.  Note  we  describe 
the  rectangle  and  trapezoid  functions  using  sums  of  unit  step  and  ramp  functions. 

h0(x,y,z,t)  =  A0(x,y,z)8(t-t0)i  (3) 
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K (*, y,z,t)  =  Al (x, y, z)u (f -  tl0 )  -  A{(x,y,z)u(t  -  tx , )  (4) 

h2 (x,y,z,t)  =  (t- 120 )A2(x,y,z)u(t-t20)-(t-t2  i)A2(x,y,z)u(t-t2  l) 
-(t-t2.2)A2(x,y,z)u(t-t22)  +  (t-t23)A2(x,y,z)u(t-t23) 


where  h0,  h{ ,  and  h2  represent  the  delta,  rectangle,  and  trapezoid  spatial  impulse 
responses  respectively  and  u(t)  is  the  unit  step  function.  The  scaling  functions 
A0(x,y,z),  A,(x,y,z),  and  A1(x,y,z)  are  constant  at  any  specific  spatial  location  and 
include  1/r  spreading,  scaling  to  account  for  the  element  size,  and  an  obliquity  factor  to 
account  for  a  soft  transducer  baffle  [8],  if  desired.  The  time  delay  t0  present  in  (3)  is 
determined  by  the  speed  of  sound  and  the  distance  from  the  element  center  to  the  field 
point.  Similarly,  the  delays  present  in  (4)  and  (5)  are  determined  from  the  speed  of  sound 
and  distances  between  the  field  point  and  the  element  edges  and  comers,  respectively. 

To  determine  the  two-way  response  from  an  element  pair  we  must  convolve  the 
appropriate  version  of  (3)-(5)  for  the  receive  element  with  the  appropriate  version  of  (3)- 
(5)  for  the  transmit  element.  While  superficial  analysis  suggests  that  nine  permutations 
are  possible,  a  more  careful  examination  reveals  that  since  the  order  of  convolution  is 
irrelevant,  some  of  these  permutations  are  redundant.  Thus,  the  two-way  response  must 
fit  one  of  the  following  six  general  expressions. 


h,r  -  h,  *hr  =  | hfi *h0  or  h0 */z,  or  h0*h2  or  /?,  *hl  or  \  *h2  or  h2  *h2 


(6) 


where  we  have  dropped  space  and  time  references  to  simplify  notation.  Note  that  while 
the  simplified  notation  of  (6)  suggests  that  in  some  cases  (such  as  h0*h0 )  the  transmit 

and  receive  responses  are  identical,  this  is  intended  only  to  state  that  the  transmit  and 
receive  responses  fit  the  same  function;  they  may  have  different  delays  and  scaling. 
Substituting  (3)-(5)  into  (6)  yields  a  set  of  six  possible  two-way  impulse  responses: 


Vfy  A)A^(  l)  to  hj)  (8) 

j= 0 


A ^0  2,j)U{t  t0  t2j  j  (9) 

j= 0 


1  1 


K*K  A0A(,XX(  0  (  !)  (f  *1  h.,j  hb.k) 


j=0 k=0 


(10) 
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i  3 

*^2  —  0  ck (t  —  tt  j  —  f2t*)  u[t  —  —  t2k^j  (11) 

/=0  *=o 

3  3  3 

\  =  (*  ~  ?2„,y  ~  f2„,*)  M(*  ~  {2 „,j  ~  *2 „,*)  (12) 

j=0k=0 

where  c;  =  {l  for  j  =  0  or  3  and  - 1  /or  j  =  1  or  2} .  With  each  of  the  six  possible 

two-way  impulse  responses  in  hand  we  can  now  consider  an  appropriate  analytical 
representation  of  the  excitation  function  (including  transmit  and  receive  electromagnetic 
transfer  functions).  While  a  number  of  possible  functions  are  attractive,  we  choose  to 
represent  emmlr(t )  using  cubic  splines  [9].  This  representation  is  attractive  because  it 
allows  arbitrary  function  shapes  while  restricting  the  form  of  the  function  to  be  no  higher 
order  than  piecewise  cubic  polynomial.  Writing  this  spline  representation  explicitly 
yields: 


emm,r(t)=  ^  (cCj  +  p/ +  yf  +  8/)(u(t  -  jd)-u(t -{j +  l)d))  (13) 

j=M„ 

where  aj ,  ft . ,  y j ,  and  Sj  are  the  spline  coefficients,  M0  and  M,  are  the  first  and  last 
spline  indices,  and  dj  is  the  spline  to  spline  interval.  We  can  now  complete  an  analytical 

expression  for  (2)  by  convolving  (13)  with  the  appropriate  version  of  (7)-(12).  While 
such  a  convolution  appears  quite  tedious,  it  can  be  readily  performed  utilizing  Laplace 
transforms  [10].  The  resulting  expression  is  a  sum  of  3rd,  4th,  5th,  6th,  or  7th  order 
polynomials  (multiplied  by  unit  step  functions),  with  the  polynomial  order  depending 
upon  the  combination  of  element  responses  employed.  We  do  not  present  the  analytical 
form  of  these  expressions  since  they  offer  little  insight  and  are  quite  lengthy.  Utilizing 
these  final  polynomial  expressions,  the  complete  two-way  response  for  a  given  transmit- 
receive  element  pair  can  be  computed  by  simply  summing  polynomials. 

Validation: 

The  proposed  algorithm  was  implemented  in  a  C  routine  called  as  a  mex  file 
within  Matlab  (The  Mathworks,  Inc.  Natick,  MA).  This  approach  allowed  us  to  readily 
generate  array  geometries  and  visualize  results  within  Matlab,  while  taking  advantage  of 
the  increased  computational  efficiency  of  compiled  C.  All  calculations  were  performed  in 
IEEE  Standard  double  precision  floating  point  arithmetic.  All  simulations  were 
performed  on  an  Apple  Powerbook  G4  computer  running  Matlab  7.0.1  under  Mac  OS 
10.3.7. 

The  validity  of  the  proposed  algorithm  was  tested  by  comparing  the  two-way 
spatial  response  of  a  single  2D  array  element  as  predicted  by  DELFI  and  FIELD  II.  For 
both  codes  the  array  element  was  modeled  using  a  single  computational  element.  Both 
codes  also  assumed  a  combined  excitation  and  transmit/receive  electromechanical 
impulse  response  equal  to  a  5.0  MHz  sine  multiplied  by  an  8  cycle  Nuttall  window  [11]. 
The  spline  representation  of  the  pulse  utilized  in  DELFI  was  generated  at  a  sampling  rate 
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of  8  times  the  assumed  center  frequency  (40  MHz).  This  sampling  rate  was  found 
empirically  to  generate  a  pulse  with  harmonics  approximately  100  dB  down  from  the 
main  signal.  FIELD  II  calculations  utilized  a  sampling  rate  5  times  higher  than  this  (200 
MHz).  While  this  rate  is  twice  the  default  for  FIELD  II  (100  MHz),  our  experience 
suggests  that  this  rate  is  certainly  within  an  appropriate  range  for  calculations  of  this  type. 
The  system  response  was  determined  in  radial  coordinates  over  an  angle  of  0°  to  90° 
(sampled  at  0.45°)  covering  a  range  from  19.925  cm  to  20.075  cm  (sampled  at  10  pm). 
For  DELFI  all  points  were  computed  at  the  time  instant  that  centered  the  sensitivity 
function  at  20  cm.  For  FIELD  II  a  complete  response  in  range,  azimuth,  and  time  was 
computed,  with  the  time  that  centered  the  sensitivity  function  at  20  cm  selected  for 
analysis. 

Typical  DELFI  and  FIELD  II  sensitivity  functions  for  a  300  x  300  pm  element 
are  shown  in  figure  2.  Both  responses  have  been  normalized  to  allow  comparison.  These 
responses  are  quite  similar  except  for  a  series  of  artifacts  located  about  the  0°  line  in  the 
FIELD  II  response.  Although  the  source  of  this  artifact  is  not  apparent,  it  may  result  from 
the  corrections  FIELD  II  employs  when  sampling  the  infinite  bandwidth  spatial  impulse 
responses  of  the  array  elements. 

Additional  simulations  were  performed  to  determine  computational  times  and  the 
similarity  of  responses  for  square  elements  ranging  in  size  from  50  to  600  microns.  All 
computational  times  were  estimated  using  the  tic  and  toe  commands  in  Matlab.  Across 
this  set  of  12  simulations  the  computation  time  for  DELFI  ranged  from  approximately 
0.21  to  0.29  seconds.  Over  the  same  set  of  conditions  FIELD  II  execution  times  ranged 
from  5.33  to  7.61  seconds.  Comparing  matched  sets  of  simulations,  FIELD  II  took  from 
18.4  to  33.6  times  as  long  for  the  computation,  with  the  average  time  of  computation 
being  27.1  times  greater  for  FIELD  II.  Note  that  these  times  do  not  include  any  array 
definitions  or  other  housekeeping  operations.  For  the  sake  of  comparison  we  also 
computed  the  temporal  response  at  the  focus  using  both  codes.  FIELD  II  computations 
were  performed  at  a  sampling  rate  of  200  MHz  and  the  DELFI  code  was  run  repeatedly  at 
different  time  points  to  generate  a  waveform  sampled  at  40  MHz.  Under  these  conditions 
the  FIELD  II  response  was  computed  in  1.2  ms  and  the  DELFI  response  was  computed  in 
7.2  ms,  giving  FIELD  II  a  factor  of  6  advantage  in  computational  time.  This  difference  is 
surprisingly  modest  given  that  DELFI  was  called  65  times  from  within  a  loop  in  Matlab, 
a  process  that  engendered  significant  overhead.  Straightforward  modifications  of  DELFI 
might  significantly  improve  its  relative  performance  when  temporal  responses  are 
required.  While  it  would  be  inappropriate  to  over-generalize  from  this  limited  set  of 
simulations,  the  existing  DELFI  code  is  computationally  attractive  for  applications  where 
spatial  responses  at  a  single  time  point  are  desired. 

The  accuracy  of  DELFI  was  quantified  by  computing  the  Full  Width  at  Half 
Maximum  (FWHM)  of  the  two-way  single  element  sensitivity  function  over  a  range  of 
element  sizes  ranging  from  50  to  600  microns  in  50  micron  steps.  The  same  metric  was 
also  computed  for  FIELD  II  simulations,  with  both  results  compared  to  the  iterative 
theoretical  method  of  [12].  McKeighen’s  iterative  method  determines  the  element’s 
FWHM  sensitivity  by  modeling  the  response  as  a  sine  squared  multiplied  by  a  cosine 
squared  obliquity  term  (to  account  for  an  assumed  soft  baffle  condition).  McKeighen’s 
method  was  implemented  using  100  iterations  at  an  assumed  center  frequency  of  5.0 
MHz.  FIELD  II  and  DELFI  used  sampling  rates  and  the  pulse  definitions  described 
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above.  The  FWHM  of  the  spatial  sensitivity  function  was  estimated  from  the  curve 
formed  by  taking  the  maximum  of  each  lateral  line,  after  computing  the  absolute  value  of 
the  response.  Results  comparing  FIELD  II,  DELFI,  and  McKeighen’s  theory  are  shown 
in  figure  3.  All  three  methods  are  in  excellent  agreement. 

Discussion: 

Computing  the  spatial  impulse  response  of  a  transmit-receive  element  pair  by 
convolving  (13)  with  (7)-(12)  analytically  offers  both  advantages  and  disadvantages 
relative  to  the  established  approach  using  discrete  time  convolution.  On  the  positive  side 
the  analytical  approach  allows  direct  computation  of  the  response  at  a  single  instant  in 
time.  In  contrast,  the  conventional  discrete  time  approach  (at  least  as  currently 
implemented)  requires  the  computation  of  a  full  temporal  response  even  if  only  a  single 
time  point  is  required.  An  additional  advantage  of  the  analytical  approach  is  that  it  yields 
an  exact  solution,  at  least  within  the  numerical  precision  of  the  computer  used  for  the 
calculation  and  within  the  limitations  of  the  far-field  approximation.  In  contrast,  the 
discrete  time  approach  utilizes  under-sampled  versions  of  the  element  responses  and  thus 
introduces  some  error.  In  addition,  the  discrete  implementation  makes  it  difficult  to 
compute  the  spatial  response  at  any  exact  instant  in  time.  A  relative  weakness  of  the 
analytical  approach  is  that  significantly  more  computation  is  required  to  determine  a 
temporal  response  at  a  single  location.  A  further  limitation  of  the  analytical  approach,  at 
least  as  implemented  here,  is  the  occasional  numerical  instability  that  results  from 
subtracting  multiple  high-order  polynomials. 

The  DELFI  code  described  here  is  still  in  an  early  stage  of  development,  with 
many  directions  apparent  for  enhancing  computational  efficiency  and  improving 
accuracy.  The  current  code  determines  the  result  of  the  convolution  of  a  single  spline 
segment  with  two  trapezoid  functions  by  summing  16  7th  order  polynomials.  This 
approach  is  computationally  tedious  and  prone  to  numerical  instability.  One  possible 
direction  for  simplification  would  begin  with  the  recognition  that  the  first  half  of  this 
response  (in  time)  requires  only  8  polynomials  for  synthesis.  Next,  by  considering  the 
response  in  negative  time  it  is  apparent  that  the  second  half  (again  in  time)  of  this 
response  also  requires  only  8  polynomials.  Together  this  realization  could  cut 
computation  time  in  half.  Additional  computational  savings  are  undoubtedly  possible. 
Algorithmic  accuracy  could  be  enhanced  by  utilizing  a  more  sophisticated  model  for  the 
element  impulse  responses.  Spline  functions  might  prove  particularly  suitable  for 
modeling  complicated  near-field  responses  [13,  14]. 

Conclusion: 

The  DELFI  code  presented  here  employs  analytical  convolution  of  cubic  spline 
functions  with  continuous  time  element  responses  to  compute  the  spatial  impulse 
response  of  ultrasound  transducers.  Comparison  with  FIELD  II  and  theoretical  methods 
show  that  the  proposed  algorithm  performs  comparably  to  existing  methods.  For  the 
computation  of  impulse  responses  at  a  single  instant  in  time,  the  proposed  algorithm  is  as 
much  as  30  times  faster  than  FIELD  II.  For  computing  temporal  responses  the  proposed 
code  is  approximately  6  times  slower  than  FIELD  II,  although  further  refinement  may 
reduce  this  mismatch.  The  proposed  algorithm,  as  implemented  by  the  DELFI  code, 
offers  an  attractive  complement  to  the  well  established  FIELD  II  program. 
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Figure  1: 

Geometries  for  determining  the  one-way  spatial  impulse  response  of  an  individual  array 
element.  In  the  left  panel  the  field  point  lies  on  a  line  through  the  element’s  center  and 
perpendicular  to  its  face.  In  the  center  panel  the  field  point  does  not  satisfy  the  first 
condition,  but  lies  on  a  plane  that  bisects  the  element  and  is  perpendicular  to  its  face.  In 
the  right  panel  the  field  point  lies  at  any  location  not  satisfying  either  of  the  first  two 
conditions. 
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Figure  2: 

The  spatial  impulse  response  of  a  300  pm  square  array  element  in  a  soft  baffle  for  an  8 
cycle  Nuttall  windowed  5.0  MHz  transmit  pulse.  The  left  panel  depicts  the  output  of  the 
proposed  algorithm  (DELFI)  while  the  right  panel  indicates  the  output  of  FIELD  II.  Both 
responses  have  been  normalized  to  allow  easy  comparison.  The  responses  are  nearly 
identical  with  the  exception  of  localized  artifacts  near  the  0°  line  of  the  FIELD  II 
response. 
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Element  Width  (urn) 


Figure  3: 

Full  Width  at  Half  Maximum  of  the  spatial  sensitivity  function  of  a  square  2D  array 
element.  All  analyses  /  simulations  were  performed  at  a  5.0  MHz  center  frequency  with 
the  soft  baffle  assumption.  All  methods  are  clearly  in  excellent  agreement. 
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